/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#define max(x,y) (((x)>(y))? (x):(y))
#define min(x,y) (((x)<(y))? (x):(y))

#include "linloglayout.h"
#include "graph.h"
#include "controller.h"

LinlogLayout::LinlogLayout(const char *name): Layout("ForceDirected", name) {
  setTimestep(1.0);
  setRunning(false);
  timestep = 0.1;
  T = 0.1; // temperature

  setPropertyFloat("Time step", 1.0);
  setPropertyFloat("Temperature", T);
  setPropertyBool("Enabled", false);
  setPropertyBool("Reset", false);
  setPropertyInt("Iterations", 100);
  setPropertyFloat("EdgeThreshold", 0.0);

  setPropertyFloat("Factor.Repulsion", 0.1);
  setPropertyFloat("Factor.Attraction", 0.1);
  setPropertyFloat("Factor.Gravitation", 0.1);

  setPropertyStringList("EdgeThresholdProperty", "Multiplicity", "edge");
  setPropertyStringList("Weight.Edge", "Default", "edge");
  setPropertyStringList("Weight.Node", "Default", "node");

  timer = createTimer(20, this);  
}

void LinlogLayout::setTimestep(float ts) {
  timestep = ts;
}

float LinlogLayout::getEnergy() {
  return totalEnergy;
}

bool LinlogLayout::isRunning() {
  return running;
}

void LinlogLayout::setRunning(bool value) {
  this->running = value;
}

float randomf2(float a, float b) {
  return (a +  b * (float) rand()/(float) RAND_MAX);
}

void LinlogLayout::updateVisible() {
  visible.clear();
  for(int i=0;i<numNodes();i++) {
    if(isNodeinLayout(i)) {
      visible.push_back(i);
    }
    else setForce(i, 0, 0);
   }
}

#define sgn(x) (((x)<0)? (-1.0):(1.0))
void LinlogLayout::setGraph(Graph *g) {
  Layout::setGraph(g);
  energy.clear();
  for(unsigned int i=0;i<graph->nodes.size();i++) {
    energy.push_back(1e5);
  }
}

float LinlogLayout::getDist(int i, int j) {
  float xi, yi, xj, yj;
  getLocation(i, xi, yi);
  getLocation(j, xj, yj);
  return sqrt((xj-xi)*(xj-xi) + (yj-yi)*(yj-yi));
}

float LinlogLayout::getDistToBarycenter(int i) {
  float xi, yi, xj, yj;
  getLocation(i, xi, yi);
  xj = barycenter[0];
  yj = barycenter[1];
  return sqrt((xj-xi)*(xj-xi) + (yj-yi)*(yj-yi));
}

void LinlogLayout::getDirection(int i, float *dir) {
  dir[0] = dir[1] = 0;
  float dir2 = 0;
  float x, y;
  getLocation(i, x, y);
  string property = getPropertyString("Weight.Edge");
  string propertyNode = getPropertyString("Weight.Node");
  float nodeweight = (graph->nodes[i]->hasProperty(propertyNode))? graph->nodes[i]->getPropertyFloat(propertyNode): 1.0;

  // Attraction Dir
  for(int k=0;k<numNeighbors(i);k++) {
    int j = getNeighbor(i,k);
    if(isEdgeinLayout(i,j)) {
      float xj, yj;
      getLocation(j, xj, yj);
      float d = getDist(i,j);
      if(d==0) continue;
      Edge * e = graph->nodes[i]->getEdge(j);
      float weight = (e->hasProperty(property))? fabs(e->getPropertyFloat(property)): 1.0;

      float tmp = weight*pow(d, attrExponent-2);
      //float tmp = weight*1.0/d;//pow(d, attrExponent-2);
      dir2 += tmp*fabs(attrExponent-1);
      dir[0] += (xj-x)*tmp;
      dir[1] += (yj-y)*tmp;
    }
  }
  
  // Repulsion Dir
  for(unsigned int k=0;k<visible.size();k++) {
    int j = visible[k];
    if(i==j) continue;
    float nodeweight2 = (graph->nodes[j]->hasProperty(propertyNode))? graph->nodes[j]->getPropertyFloat(propertyNode): 1.0;
    float xj, yj;
    getLocation(j, xj, yj);
    float d = getDist(i,j);
    float tmp = repuFactor * nodeweight * nodeweight2 * pow(d, repuExponent-2);
    dir2 += tmp* fabs(repuExponent-1);
    dir[0] -= (xj-x)*tmp;
    dir[1] -= (yj-y)*tmp;
  }

  // Gravitation Dir 
  double d = getDistToBarycenter(i);
  double tmp = gravFactor * repuFactor * nodeweight * pow(d, attrExponent-2);
  dir[0] += (barycenter[0] - x)*tmp;
  dir[1] += (barycenter[1] - y)*tmp;
  dir2 += tmp * fabs(attrExponent-1);

  // compute average Euclidean distance to other nodes
  float avgDist = 0.0;
  for(unsigned int k=0;k<visible.size();k++) {
    int j = visible[k];
    avgDist += getDist(i,j);
  }
  avgDist /= (float) ((int) visible.size()-1);
  if (dir2 != 0.0) {
    // normalize force vector with second derivation of energy
    dir[0] /= dir2;
    dir[1] /= dir2;
    // ensure that the length of dir is not greater
    // than average Euclidean distance to other nodes
    double length = sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
    if (avgDist > 0.0 && length > avgDist) {
      length /= avgDist;
      dir[0] /= length;
      dir[1] /= length;
    }
  } else {
    dir[0] = dir[1] = 0.0;
  }
}

float LinlogLayout::getEnergy(int i){
	string edgeFilterProperty = PROPERTY(getPropertyString("Edge Filter"));
  float energy = 0;
  string property = getPropertyString("Weight.Edge");
  string propertyNode = getPropertyString("Weight.Node");
  float nodeweight = (graph->nodes[i]->hasProperty(propertyNode))? graph->nodes[i]->getPropertyFloat(propertyNode): 1.0;
  // Attraction
  float totalAttraction = 0;
  float totalRepulsion = 0;
  for(int k=0;k<numNeighbors(i);k++) {
		bool useEdge = graph->getAsBoolean(graph->nodes[i]->edges[k], edgeFilterProperty);
    int j = getNeighbor(i,k);
    if(isEdgeinLayout(i,j) && useEdge) {
      float d = getDist(i,j);
      if(d==0.0) continue;
      Edge * e = graph->nodes[i]->getEdge(j);
      float weight = (e->hasProperty(property))? fabs(e->getPropertyFloat(property)): 1.0;
      float nodeEnergy;
      if (attrExponent == 0.0) {
	nodeEnergy = weight*log(d);
      } else {
	nodeEnergy = weight*pow(d, attrExponent) / attrExponent;
      }
      //if(isnan(nodeEnergy)) printf("Attr.Energy[%d-%d] = %f d = %f\n", i,j, nodeEnergy, d);
      energy += nodeEnergy;
      totalAttraction += nodeEnergy;
    }
  }
  
  // Repulsion
  for(unsigned int k=0;k<visible.size();k++) {
    int j = visible[k];    
    if(i==j) continue;
    float nodeweight2 = (graph->nodes[j]->hasProperty(propertyNode))? graph->nodes[j]->getPropertyFloat(propertyNode): 1.0;
    float d = getDist(i,j);
    float nodeEnergy;
    if(d==0.0) continue;
    if (repuExponent == 0.0) {
      nodeEnergy = repuFactor*nodeweight*nodeweight2*log(d);
    } else {
      nodeEnergy = repuFactor*nodeweight*nodeweight2*pow(d, repuExponent) / repuExponent;
    }
    //if(isnan(nodeEnergy)) printf("Rep.Energy[%d-%d] = %f d = %f\n", i,j, nodeEnergy, d);
    totalRepulsion += nodeEnergy;
    energy -= nodeEnergy;
  }


  // Gravitation
  double d = getDistToBarycenter(i);
  if (attrExponent == 0.0) {
    energy += gravFactor * nodeweight * log(d);
  } else {
    energy += gravFactor * nodeweight * pow(d, attrExponent) / attrExponent;
  }

  //printf("energy[%d] = %f\n", i, energy);
  return energy;
}

void LinlogLayout::computeBaryCenter() {
  barycenter[0] = barycenter[1] = 0.0;
  double weightSum = 0.0;
  float x, y;

  for(unsigned int ki=0;ki<visible.size();ki++) {
    int i = visible[ki];
    string propertyNode = getPropertyString("Weight.Node");
    float nodeweight = (graph->nodes[i]->hasProperty(propertyNode))? graph->nodes[i]->getPropertyFloat(propertyNode): 1.0;

    
    weightSum += nodeweight;
    getLocation(i, x, y);
    barycenter[0]+=nodeweight*x;
    barycenter[1]+=nodeweight*y;
  }
  if (weightSum > 0.0) {
    barycenter[0]/=weightSum;
    barycenter[1]/=weightSum;
  }
}

int step;
void LinlogLayout::computeIteration() {
  if(!running || timestep<=0) return;
  //float edgeThresh = getPropertyFloat("EdgeThreshold",0);
  int nrIterations = getPropertyInt("Iterations");

  float finalAttrExponent = 1.0; // LIN
  float finalRepuExponent = 0.0; // LOG

  float totalEnergy = 0;
  float minx, maxx, miny, maxy;
  minx = miny = 9999.99;
  maxx = maxy = -9999.99;

  float bestDir[2];
  computeBaryCenter();


  // Cooling off
  attrExponent = finalAttrExponent;
  repuExponent = finalRepuExponent;
  if(finalRepuExponent<1.0) {
    if (step <= 0.6*nrIterations) {
      // use energy model with few local minima 
      attrExponent += 1.1 * (1.0 - finalRepuExponent);
      repuExponent += 0.9 * (1.0 - finalRepuExponent);
    } else if (step <= 0.9*nrIterations) {
      // gradually move to final energy model
      
      attrExponent += 1.1 * (1.0 - finalRepuExponent)
	* (0.9 - ((double)step)/nrIterations) / 0.3;
      repuExponent += 0.9 * (1.0 - finalRepuExponent)
	* (0.9 - ((double)step)/nrIterations) / 0.3;
      
    }
  }
  //attrExponent = finalAttrExponent;
  //repuExponent = finalRepuExponent;
  printf("Exponents %f %f\n", attrExponent, repuExponent);
  step++;

  for(unsigned int ki=0;ki<visible.size();ki++) {    
    if(ki%50==0) { printf("."); fflush(stdout);}
    int i = visible[ki];
    if(graph->nodes[i]->getPropertyFloat("Degree")==0) continue;
    float oldEnergy = getEnergy(i);

    // compute direction of the move of the node
    getDirection(i, bestDir);
    float x, y;
    float x0, y0;
    getLocation(i, x, y);
    x0 = x;
    y0 = y;   

    float bestEnergy = oldEnergy;
    float bestMultiple = 0;
    bestDir[0]/=32.0f;
    bestDir[1]/=32.0f;
    for (int multiple = 32;  multiple >= 1 && (bestMultiple==0 || bestMultiple/2==multiple); multiple /= 2) {
      x = x0 + bestDir[0] *(float) multiple;
      y = y0 + bestDir[1] *(float) multiple;
      setLocation(i, x, y);
      float curEnergy = getEnergy(i);
      if (curEnergy < bestEnergy) {
	bestEnergy = curEnergy;
	bestMultiple = multiple;
      }
    }
    for (int multiple = 64; multiple <= 128 && bestMultiple == multiple/2;  multiple *= 2) {
        x = x0 + bestDir[0] *(float) multiple;
	y = y0 + bestDir[1] *(float) multiple;
	setLocation(i, x, y);
 	float curEnergy = getEnergy(i);
	if (curEnergy < bestEnergy) {
	  bestEnergy = curEnergy;
	  bestMultiple = multiple;
	}
    }

    x = x0 + timestep*bestDir[0] * bestMultiple;
    y = y0 + timestep*bestDir[1] * bestMultiple;
    setLocation(i, x, y);
    totalEnergy += bestEnergy; 
    graph->nodes[i]->setPropertyFloat("Energy", bestEnergy);
  }
  T = T*0.99; // cooling off
  print("\n[%d] E = %lf\n", step, totalEnergy);
	updated();
}

void LinlogLayout::initEnergyFactors() {
  attrExponent = 1.0;
  repuExponent = 0.0;
  string property = getPropertyString("Weight.Edge");
  string propertyNode = getPropertyString("Weight.Node");

  for(unsigned int ki=0;ki<visible.size();ki++) {
    int i = visible[ki];
    graph->nodes[i]->setPropertyFloat("Energy", 1e10);
  }

  double attrSum = 0.0;
  for(unsigned int ki=0;ki<visible.size();ki++) {
    int i = visible[ki];
    if(graph->nodes[i]->getPropertyFloat("Degree")==0) continue;
    for(int k=0;k<numNeighbors(i);k++) {
      int j = getNeighbor(i,k);
      if(isEdgeinLayout(i,j)) {
	Edge * e = graph->nodes[i]->getEdge(j);
	float weight = (e->hasProperty(property))? fabs(e->getPropertyFloat(property)): 1.0;
	attrSum += weight;
      }
    }
  }
  double repuSum = 0.0;
  for(unsigned int ki=0;ki<visible.size();ki++) {
    int i = visible[ki];
    if(graph->nodes[i]->getPropertyFloat("Degree")==0) continue;
   float nodeweight = (graph->nodes[i]->hasProperty(propertyNode))? graph->nodes[i]->getPropertyFloat(propertyNode): 1.0;
    repuSum += nodeweight;
  }

  if (repuSum > 0.0 && attrSum > 0.0) {
    double density = attrSum / repuSum / repuSum; 
    repuFactor = density * pow(repuSum, 0.5*(attrExponent-repuExponent));
    gravFactor = density * repuSum * pow(gravFactor, attrExponent-repuExponent); 
  } else {
    repuFactor = 1.0;
  }
  printf("Attr.sum = %f repuSum = %f\n", attrSum, repuSum);
  printf("Energy Factors = %f %f\n", repuFactor, gravFactor);
}

int LinlogLayout::numNodes() {
  return graph->nodes.size();
}

inline bool LinlogLayout::isNodeinLayout(int i) {
  return graph->nodes[i]->isVisible();
}
bool LinlogLayout::isEdgeinLayout(int i, int j) {
	return isNodeinLayout(i) && isNodeinLayout(j);
}
bool LinlogLayout::isNodeUpdatable(int i) {
  return !graph->nodes[i]->isSelected(); 
}
int LinlogLayout::numNeighbors(int i) {
  return graph->nodes[i]->edges.size();
}
int LinlogLayout::getNeighbor(int i, int k) {
  return graph->nodes[i]->edges[k]->getNode1()->getIndex();
}

void LinlogLayout::reset() {
  for(unsigned int i=0;i<graph->nodes.size();i++) {    
    float x,y;
    x = randomf2(-1,1);
    y = randomf2(-1,1);
    setLocation(i,x,y);
    setVelocity(i,0,0);
  }
}

void LinlogLayout::timerEvent() {
    updateVisible();
    computeIteration();
    controller->updateRender();
}

void LinlogLayout::compute() {
  int numIterations = getPropertyInt("Iterations");
  timestep = getPropertyFloat("Time step");
  step = 0;
  updateVisible();
  initEnergyFactors();
  running = true;
  for(int i=0;i<numIterations;i++) {
    listener->setProgress((float) i/(float) numIterations);
    updateVisible();
    computeIteration();
  }
  running = false;
}

void LinlogLayout::update(string property) {
	Layout::update(property);
  timestep = getPropertyFloat("Time step");
  float temp = getPropertyFloat("Temperature");
  T = temp;
  if(getPropertyBool("Reset")) {
    reset();
  }
  bool enabled = getPropertyBool("Enabled");
  if(enabled) {
    if(!running) {
      step = 0;
      updateVisible();
      initEnergyFactors();
      timer->start();
      running = true;
    }
  } else {
    if(running) {
      timer->stop();
      running = false;
    }
  }
}

Layout *createLinlogLayout(string name) {
  LinlogLayout *layout = new LinlogLayout(name.c_str());
  return layout;
}
